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Abstract 



To calculate dilepton rates in a Monte Carlo simulation of ultrarelativistic 
heavy ion collisions, one usually scales the number of similar QCD processes by 
a ratio of the corresponding differential probabilities. We derive the formula 
for such a ratio especially for dilepton bremsstrahlung processes. We also 
discuss the non-triviality of including higher order corrections to direct Drell- 
Yan process. The resultant mass spectra from our Monte Carlo simulation are 
consistent with the semi-analytical calculation using dilepton fragmentation 
functions. 
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I. INTRODUCTION 



In a previous paper [0, we investigated dilepton production associated with minijet 
final state radiation in heavy ion colhsions at colhder energies, using dilepton fragmentation 
functions which can be evaluated perturbatively. The dilepton pairs from the fragmentation 
of minijets are found to be comparable to direct Drell-Yan (DY) for small invariant mass M ~ 
1-2 GeV/c^ at the highest energy of the Brookhaven Relativistic Heavy Ion Collider (RHIC). 
At the CERN Large Hadron Collider (LHC) energies, the associated dilepton production 
becomes dominant over a relative large range of the invariant mass. 

Due to the relatively large invariant dilepton mass, M ^ A, the radiative corrections 
are calculable in pQCD up to all orders in the leading logarithm approximation. Collinear 
approximation is also used in convoluting the obtained dilepton fragmentation functions with 
minijet cross sections to compute the radiative contributions to dilepton production. Since 
there exist Monte Carlo simulations of QCD cascading pj-Q which can take into account 
many other effects, like multiple ladder structure, it is important to check our semi-analytical 
approach with realistic Monte Carlo simulations. In this way, we can address the validity of 
the approximations we made in the fragmentation function approach 

To directly simulate dilepton production in a Monte Carlo event generator is rather 
difficult due to the small QED coupling constant as compared to that of QCD. To overcome 
this difficulty, one can multiply the number of specific QCD processes, which resemble those 
of dilepton production by an appropriate ratio of the corresponding differential probabilities, 
as has been tried by Geiger and Kapusta However, the problem is more complicated than 
one might first think. For radiative dilepton production, one has to take into account the fact 
that the corresponding radiated quarks and gluons in QCD can have further bremsstrahlung 
which is different from the QED case. The additional bremsstrahlung gives rise to an extra 
Sudakov form factor and one must include it in the differential ratio to give the correct 
dilepton emission. We will derive a formula for the differential ratio and demonstrate that 
the resultant simulation is consistent with our previous semi-analytical calculation in Ref . [|1| . 
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The problem of simulating direct DY process among QCD processes lies in how to take 
into account higher order, 0{a'^as), corrections. We will discuss when an overall K factor, 
as used in most Monte Carlo simulations of QCD hard processes in hadronic and nuclear 
collisions, is sufficient enough to simulate the QCD contributions. We also consider how 
possible double counting can be avoided. 

The remainder of the paper is organized as follows. In the next Section, we will dis- 
cuss how to calculate the differential ratio to obtain dilepton emission from a Monte Carlo 
simulation of the corresponding pQCD processes. In Section III, we will discuss how to 
simulate direct DY processes, especially how to take into account higher order corrections. 
The results of our simulation will be compared to semi-analytical calculations in Section IV, 
and finally, a summary with some discussions is given in Section V. 

II. SIMULATION OF DILEPTON BREMSSTRAHLUNG 

If a Monte Carlo generator does not have the QED processes built in, one can calculate 
the absolute cross sections of the QED processes by scaling the differential cross sections of 
certain types of QCD processes. The scaling factor however depends on both the QED and 
QCD processes. For dilepton production through bremsstrhlung, one would naively think 
that a scaling factor between virtual photon and virtual gluon radiation processes from a 
quark line is enough. However, the probability to find a virtual gluon with fixed invariant 
mass depends on the probability that the gluon does not have any further radiations to 
degrade its virtuality. One therefore should use the scaling factor between process (a) and 
processes (b) and (c) in Fig. |l[ 

The Monte Carlo simulation of QCD cascading is carried out by giving for each vertex of 
the radiation tree, such as those in Fig. |l|(b), a norma/ueo? probability distribution. Given 
the maximum virtuality Q^ax of ^ particular process, the normalized probability for the 
off-shell parton a with < Q^ax branch into partons b and c of light-cone momentum 
fractions z and 1 — z is [0,^ , 
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The relative transverse momentum between the radiated partons h and c is given by, 

g^ = .(l-z)(g^-|-^). (2) 

Note that the variable z{l — z)q^ in the strong coupling constant in Eq. |] is approximately 
In a Monte Carlo simulation, the time- like branching is usually terminated at q^ < /Xq, 
where the physics of nonperturbative hadronization sets in. By requiring g^, ^ /^o ^^"^ 
relative transverse momentum qx to be real, one defines the kinematically allowed region of 
the phase space as 

4/^0 < (f< <5max; 

e(g) <z<l- e(g), e(g) = \{l- ^l-A^^l/q^). (3) 

If the virtuality of one of the radiated partons is fixed to q^ = M^, the above region of the 
phase space is modified to. 



ei{q,M) < z<e2{q,M), (4) 

where, 



ei,2(g,M) = i 



\ 



Q J Q 



(5) 



In Eq. |T], the Sudakov form factor Saiq"^) is defined as 



so that Sa{Q\^.^^) / Sa{q^) is the probability for parton a not to have any branching between 
Q^^jjx ^"^^ ■ The contribution of QED processes to the Sudakov form factor is negligible. 

Based on the above probability distributions, we can write down the differential proba- 
bility for a quark to radiate either a qq pair or two gluons via an intermediate virtual gluon 
with invariant mass M, as shown in Figs. |l|(a) and (b). 



max/ 



Notice the variables in the second set of Sudakov form factors. Since is the actual 
virtuality of the quark line preceding the gluon radiation and the daughter quark line has 
at least virtuality of the maximum value of the gluon virtuality, M^, is then {q — /io)^- 
Due to the same reason, the lower limit of the integration over is (M + /io)^- 

Similarly, the differential branching probability for the dilepton production via diagram 
(a) in Fig. [| is. 



dM^ M2 7(M+w)2 AifeM) "-^^"^ '27r Sg{q^) 

/ dziPy^e+i-{zi) — , (8) 
JO Ztt 

where again we neglect the QED contribution to the Sudakov form factor, and the integra- 
tion over z and Zi can be carried out analytically using Eq. ^ Notice that the differential 
probability for the dilepton bremsstrahlung has one Sudakov form factor less than the cor- 
responding QCD processes. To obtain radiative dilepton production [Fig. |l|(a)] by scaling 
the corresponding QCD processes [Figs. |l|(b) and (c)], one simply multiplies the number of 
virtual gluons from a quark line by a ratio, 

7^(M^«J.^/^, (9) 



for given and Q^ax- 

To demonstrate the effects of the Sudakov form factor on scaling QCD processes of the 
Monte Carlo simulations, we plot in Fig. |^ the ratio Tli^M"^ , Q'^^^^) (solid) as a function of 
Qmax at fixed M = 2 GeV/c^ together with the result (dashed) obtained when Sudakov form 
factors are set to unity. We also show the value of 7lo{M^, Q^ax) (dot-dashed) in which both 
the Sudakov form factors and z{l — z)q'^ dependence of the coupling constant are neglected. 
In this case the ratio is reduced to 
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which is independent of Qmax, and where 

/•l-e{M) 
Je{M) 

= 6 ln[l/e(M) - 1] + 9[e(M) - i]. (11) 

Since the Sudakov form factor takes into account additional branchings preceding the 
chosen vertex, it should suppress the probability distribution of the splitting g — > qq, gg 
at the given M. As we see in Fig. H, the ratio TZ{M^, Q^ax) is therefore enhanced relative 
to both the case when Sa = I and to 7?.o(M^, Q^ax) the inclusion of Sudakov form 
factors. The enhancement increases with Qmax as expected due to the increasing branching 
probability. Similar to the dilepton fragmentation functions, the ratio is very sensitive to 
the scale Qmax- We will discuss in Sec. ^ how we choose Qmax which is consistent with the 
Monte Carlo simulation of QCD cascading. When Sudakov form factors are set to unity, the 
dependence of the ratio on Qmax only comes from the z and dependence of the running 
strong coupling constant. If both z{l — z)q'^ and Zi{l — zi)M'^ in Eq. |^ are replaced by M^, 
the ratio TZo becomes larger and is independent of Qmax as shown in Fig. |[ 

III. SIMULATION OF DIRECT DY PROCESSES 

It is relatively easier to simulate the lowest order process of direct Drell-Yan by multi- 
plying the differential cross section of qq —>■ qiCji, gg by the ratio 

7? daqq^DY . s 

l^i=l '^^qq^qiqi "r '^^qq-'gg 

It is more subtle, however, to include QCD corrections. These higher order corrections which 
give rise to a so-called i^'-factor have been studied extensively . The problem here is how 
to include this i^'-factor in scaling QCD processes to obtain the DY cross sections. 

The first order correction to DY process in QCD comes from the "annihilation" qq — >■ 
(yf + DY, the "Compton" q+g g+DY process and the virtual corrections ^j. Like in deeply 
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inelastic lepton nucleoli scatterings (DIS), there are infrared singular and finite contributions 
from these corrections. The infrared singular and part of the finite terms can be absorbed 
into the quark and antiquark distribution functions which are defined in DIS processes and 
should be evaluated at the scale of according to Altarelli-Parisi evolution equations [0 . 
What are left over are finite and scheme-independent contributions from the higher order 
corrections. One can find detailed discussions in, e.g., Ref. [§. What we want to point out 
here is that the dominant contribution to the i^-factor of about 2 in the (pT-integrated) 
mass spectrum of the direct DY is from the virtual corrections. The contributions from real 
corrections which depend on both the quark and gluon distribution functions are relatively 
very small. Therefore, as a first approximation, we can include higher order corrections to 
direct DY process in our simulation by multiplying the QCD cross sections of qq — > qiqi,gg 
by an effective K = 2 factor. This K factor in principle now includes both real and virtual 
corrections. The quark distribution functions in the cross section should be evolved and 
evaluated at scale M^. 

In the Monte Carlo simulations, one could also include the real QCD-corrections explicitly 
to DY process by counting the number of similar QCD processes, qq ^ gg, q + g — ^ q + g and 
scaling them by some calculable ratio as has been done in Ref. 0. However, one still has 
to include the virtual corrections which can be characterized as an effective multiplicative 
factor, but which now differs from the normal overall DY fC-factor. This is exactly the 
problem one has to face if one wants to simulate the pt distribution of DY dilepton pairs, 
whose large pt tail mainly comes from real QCD-corrections. The lowest order DY process 
only contributes to the small px part of the spectrum by including the intrinsic Pt of quarks 
and anti-quarks. The virtual corrections to the lowest order DY can be taken into account 
by using an effective 'fT'-factor. However, one must be very careful not to include the real 
corrections in this effective 'i^^'-factor. 

In the Monte Carlo simulations of QCD processes in hadronic collisions, one usually also 
uses an effective K-factor to take into account higher order corrections [Q. However, this 
i^-factor should not be included when one scales the number of qq gg and q + g ^ q + g 
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processes by some ratio to calculate the real QCD-corrections to DY process. Otherwise, 
double-counting may occur. 



IV. NUMERICAL RESULTS 

To perform the Monte Carlo simulations, we use PYTHIA [^] subroutines for QCD hard 
scatterings and the associated bremsstrahlungs as adapted in HIJING model 0. HIJING 
is a Monte Carlo model developed for parton and particle production in high energy pp, pA 
and AA collisions. In this model, multiple minijet production at A^A^ level is calculated in 
the eikonal formalism fl^. As in many other models which attempt to merge low and high 
Pt dynamics, a pt cutoff scale po has to be introduced, which will limit the invariant masses 
of produced dileptons in our simulation. For nuclear interactions, binary approximation is 
assumed for independent hard scatterings. Jet quenching due to final state interaction of 
produced partons with string-like soft mean-field was also included in the original HIJING 
model [|1^]. We switch off these final state interactions to simplify our study here so that 
we can compare the numerical results with semi-analytical calculations. The initial parton 



distribution functions of a nucleon is taken to be Duke-Owens parameterization |[12| set 1. 



Nuclear shadowing and its scale dependence are also taken into account as in Ref. [|13 
Impact parameter dependence of the nuclear parton distributions is modeled in according to 
Refs. PJl^. However, at ^/s = 200 AGeV, the shadowing effects on the associated dilepton 



production are small as seen in Ref. IQ. 

During the final state radiation, we count the number of virtual gluons with given in- 
variant mass M which are radiated from quark lines. The maximum virtuality Qmax of the 
quark should be the invariant mass of its parent parton minus the minimum virtuality /iq of 
its sister parton. If there is no bremsstrahlung prior to this branching vertex, Q^^x should 
be related to the transverse momentum transfer pt of the corresponding hard scattering. 
In order to conserve both energy and momentum, the two produced partons from a hard 
scattering are combined together in PYTHIA to initiate final state radiation. The total 
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virtuality of the two parton system is chosen to be 2pT- Since both of the partons must 
have at least a virtuahty of fiQ, the maximum virtuahty of the selected quark immediately 
after the hard scattering should be Qmax = 2pT — yUo- With given and Q^ax; then can 
calculate the number of dileptons produced from the final state radiation by multiplying 
the number of these radiated virtual gluons with the ratio 7?.(M^, Q^ax)- Shown by the 
solid histogram in Fig. ^, is the invariant mass distribution of the radiated dileptons thus 
obtained for central Au + Au collisions at RHIC energy. In the simulation, the parton shower 
is terminated whenever a minimum virtuality /io = 0.5 GeV is reached. Then the parton 
is put on shell and considered real. So the minimum invariant mass of our selected virtual 
gluons, thus also of the dileptons, is 2/io = 1 GeV, according to Eq. |^. 

We also plot in Fig. ^ our semi-analytical calculation of the radiative dilepton production 
(solid curve) which agrees quite well with the simulated result. The small differences both in 
the total number and the slope of the distribution could come from several simplifications we 
made in our semi-analytical approach. As stated in Ref . |]1| , we did not fully take into account 
the kinematic restrictions (Eqs. |^, ^ at every stage of the radiation tree in the calculation of 
the dilepton fragmentation functions. The variable in the strong coupling constant is taken 
to be instead of the relative transverse momentum ^ z{l — z)q'^. The fragmentation 
function approach has only one branching tree corresponding to a simple ladder structure, 
whereas the Monte Carlo simulation takes into account all possible branching trees, thereby 
enhancing the small M dilepton production. In order to be as consistent as possible in both 
calculations in Fig. |], we have chosen the same scales Qmax = 2pT and A = 0.4 GeV in the 
dilepton fragmentation functions as have been used in the Monte Carlo simulation 

To simulate the lowest order direct DY process, we simply count the number of similar 
QCD subprocesses, qq qq, gg at fixed s = M^. We then multiply the number by the ratio 
TZby to obtain the number of direct DY dileptons, which is shown as the dashed histogram in 
Fig. ^ We also compare the result to the parton model calculation (dashed curve) with the 
same set of parton distribution functions as used in the simulation. Higher order corrections 
are included by multiplying a K = 2 factor in both the simulation and analytical calculation. 
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In terms of pt and rapidities yi^2 of minijets, the invariant mass of the dilepton is 

M2 = 2p|[l + cosh(|/i-|/2)]. (13) 

Since we have a cutoff po = 2 GeV/c for pt, the lower hmit of the dilepton mass from the 
Monte Carlo simulation is then M < 4 GeV/c^. The analytical calculation can go as low as 
the initial scale of the structure functions, Qo = 2 GeV. 

V. CONCLUSIONS 

In this paper, we have studied minijet-associated dilepton production in ultra-relativistic 
nuclear collisions through Monte Carlo simulations to check our previous semi-analytical 
calculation through the fragmentation function approach. We derived a formula for 
the differential ratio by which we multiply the number of similar QCD processes to obtain 
dilepton production from the Monte Carlo simulation of QCD cascading. Using this ratio, 
we found that our semi-analytical calculation is consistent with Monte Carlo simulations. 
The difference between the two due to some simplifications we made in the fragmentation 
function approach is small. 

Most importantly, we found that Sudakov form factors which were not included in the 
ratio in Ref. |^ are essential for us to give the right results. If neglected, the resultant 
dilepton rate from final state radiation would differ from our early semi-analytical calculation 
by orders of magnitude. Due to the same reason, the differential ratio is quite sensitive to the 
maximum virtuality Qmax of the branching processes, similar to the fragmentation function 
approach. One therefore has to choose its value to be consistent with what is used in the 
Monte Carlo simulation of QCD cascading. 

We also simulated the direct Drell-Yan processes of dilepton production and compared 
it to semi-analytical calculation in the parton model. We pointed out the complication 
in including higher order corrections in the Monte Carlo simulations and the possibility of 
double counting. We believe this is especially important when one wants to simulate dilepton 
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production through the final state parton rescatterings []5|,[l6| in a dense partonic system hke 
a quark-gluon plasma. Unlike in hadronic scatterings where infrared singularities due to 
real and virtual corrections can be absorbed into the definition of QCD evolved parton 
distribution functions, the screening mass due to resummation of hot thermal loops ||T^ 
naturally regulates the infrared divergences. However, one still has important contributions 



from both real and virtual corrections In order to take into account these corrections 
in a Monte Carlo simulation, one needs to analyze the higher order calculation in finite 
temperature QCD in detail. 
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FIGURES 

FIG. 1. Illustration of (a) dilepton, (b) quark-antiquark pair, and (c) two-gluon emission from 
a quark line. The dashed lines present the preceding radiation or scattering processes which 
kinematically determines the maximum value, Qmax) of the quark virtuality . 

FIG. 2. The ratio 7^(M^,(5max) between the probability of g — > t^l~ + q and q qiQi + q, 
gg + q processes as functions of Qmax at fixed M = 2 GeV/c^, with (solid) and without (dashed) 
the Sudakov form factors. TZq (dot-dashed) is obtained with both Sudakov form factors and the 
z{l — z)q^ dependence of neglected. A factor is divided out. 

FIG. 3. Mass spectrum of the minijet-associated (solid histogram) and direct DY (dashed 
histogram) dileptons from the Monte Carlo simulation and our direct calculation (solid and dashed 
curves, respectively) (with Qmax = 2^^, A = 0.4 GeV and jiQ = 0.5 GeV) in central Au + Au 
colUsions at = 200 AGeV. 
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